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A  technique  for  imaging  large  areas  of  the  ocean  floor  by  inverting  reverberation  data  is 
presented.  Acoustic  returns  measured  on  a  horizontally  towed  line  array  have  right-left 
ambiguity  about  the  array’s  axis  and  decreasing  cross-range  resolution  for  off-broadside  beams 
and  more  distant  scattering  sites.  By  optimal  use  of  data  taken  with  differing  array  locations  and 
orientations,  right-left  ambiguity  is  eliminated  and  resolution  is  maximized.  This  is 
accomplished  via  a  global  inversion  using  the  array’s  known  resolving  properties.  The  output  of 
the  inversion  can  be  optimally  resolved  reverberation,  scattering  coefficient,  or  physical 
properties  of  the  seafloor,  depending  upon  the  formulation.  Simulations  are  presented  to  show 
the  general  superiority  of  this  approach  to  data  averaging. 

PACS  numbers:  43.30.Gv,  43.30.Pc 


INTRODUCTION 

Our  goal  is  to  image  the  ocean  basin  using  acoustic 
reverberation  data.  When  operating  in  bottom-limited  en¬ 
vironments,  active  sonar  systems  measure  discrete  back- 
scatter  returns  from  the  seafloor.  The  nature  of  these  re¬ 
turns  depends  upon  transmission  loss  as  well  as  the 
morphology  and  geo-acoustic  properties  of  the  scattering 
site.  By  two-way  travel-time  analysis,  the  range  of  partic¬ 
ular  scattering  sites  can  be  deduced  for  respective  returns. 
Since  low-frequency  active  systems  typically  measure 
acoustic  returns  with  a  horizontally  towed  line  array,  the 
azimuth  of  scattering  sites  can  be  determined  by  beam¬ 
forming.  Beamformed  line-array  data,  however,  has  (a)  an 
inherent  right-left  ambiguity  about  the  array’s  axis,  (b) 
decreasing  angular  resolution  for  off-broadside  beams,  and 
(c)  decreasing  cross-range  resolution  for  more  distant  scat¬ 
tering  sites.  These  ambiguities  make  it  difficult  to  properly 
locate  returns  in  complex  ocean  environments  where  scat¬ 
tering  sites  are  broadly  distributed  in  range  and  azimuth.1 

For  future  reference,  we  define  an  “observation”  as  an 
active  insonification  and  subsequent  towed-array  measure¬ 
ment  of  reverberation  at  a  single  location  and  orientation. 
(We  sometimes  use  “observation”  to  refer  to  the  array 
location  in  such  a  measurement.)  We  define  a  “reverbera¬ 
tion  map”  as  the  magnitude  of  acoustic  returns  from  a 
single  observation  charted  to  the  horizontal  location  of 
their  respective  scattering  sites.  A  “scattering  strength 
map”  is  a  reverberation  map  normalized  by  source 
strength,  scattering  site  area,  and  transmission  loss. 

Averaging  together  scattering  strength  maps  has 
proven  to  be  successful  in  reducing  right-left  ambiguity, 
given  a  sufficiently  diverse  observation  geometry.2,3  How¬ 
ever,  this  technique  has  only  been  applied  to  situations 
where  isolated  scattering  features,  such  as  seamounts  or 
continental  margins,  disturb  sound  channels  that  otherwise 
have  excess  depth.  It  does  not  work  as  well  for  broadly 
distributed  backscattering  in  complex  bottom-limited 
environments.4 

Right-left  angular  ambiguity  in  ambient  noise  data  has 
been  resolved  by  an  iterative  optimization  procedure.5  This 


can  be  further  applied  to  resolve  angular  ambiguity  in 
long-range  monostatic  and  bistatic  reverberation  data  by 
collapsing  towed-array  observation  locations  to  a  single 
point  and  performing  angular  inversions  at  independent 
range  steps.6  However,  this  approach  provides  a  resolution 
no  finer  than  the  maximum  separation  between  observa¬ 
tions.  Given  spatially  distributed  observations  in  complex 
bottom-limited  environments,  higher  resolution  is  often 
necessary  for  comparisons  with  seafloor  morphology. 

We  propose  a  more  general  imaging  method  based  on 
a  global  inversion  of  spatially  distributed  observations.  In 
this  initial  formulation  we  assume  that  sound  scattered 
from  a  particular  region  is  independent  of  incident  and 
observation  angle,  in  both  the  horizontal  and  vertical.  This 
approximation  is  good  for  monostatic  and  bistatic  rever¬ 
beration  experiments  in  the  deep  ocean,  where  propagation 
angles  only  vary  over  a  narrow  vertical  width,  and  the 
azimuth  of  distant  scatterers  remains  relatively  constant 
from  observation  to  observation.  (It  can  also  be  used  to 
make  omnidirectional  estimates,  given  a  broad  distribution 
of  observations,  as  is  shown  in  a  preliminary  article.7) 

By  simulation,  we  demonstrate  the  method  for  a 
monostatic  observation  geometry  and  operational  parame¬ 
ters  used  in  bottom  reverberation  experiments  sponsored 
by  the  Office  of  Naval  Research  Special  Research  Program 
in  1991-1993.  We  first  create  a  synthetic  ocean  basin  with 
a  scattering  coefficient  representation,  and  generate  syn¬ 
thetic  reverberation  maps  for  the  given  observation  geom¬ 
etry.  We  then  attempt  to  estimate  the  true  scatter  coeffi¬ 
cients  from  the  reverberation  maps  by  (a)  a  linear  average, 
(b)  a  dB  average,  and  (c)  a  global  inversion.  Finally,  we 
test  the  accuracy  of  each  method  by  a  statistical  compari¬ 
son  of  the  estimated  and  true  scatter  coefficients.  Since 
some  methods  may  estimate  high  scatter  coefficients  better, 
we  make  the  above  comparison  separately  for  different 
magnitude  regimes. 

I.  REVERBERATION  MAPS 

To  compare  reverberation  measured  with  differing  ar¬ 
ray  positions  and  orientations,  returns  from  each  observa- 
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FIG.  1.  The  image  chosen  to  represent  the  ocean  basin.  Increasing  negative  values  indicate  increasing  depth  as  measured  from  the  ocean’s  upper  surface. 
The  data  are  extracted  and  subsampled  from  the  SRP  Geophysical  Survey9  for  a  244  x  244  km2  region.  Scattering  strength,  2[/n,/i]=  10  log(<r[/i,m]),  is 
chosen  to  be  equal  to  the  numerical  value  of  depth,  in  hundreds  of  meters,  for  respective  east  versus  north  grid  point  coordinates  m=  1,2,3,...,#*  and 
n= 1,2,3,. ..,#3,,  where  Nx=Ny=  163. 


tion  must  be  mapped  to  absolute  spatial  coordinates.  For 
real  data,  temporal  returns  measured  on  an  array  of  hy¬ 
drophones  can  be  accurately  converted  to  the  range  and 
azimuth  of  their  respective  scattering  sites.  Beamforming 
yields  the  necessary  azimuthal  dependence,  but  with  vary¬ 
ing  resolution  and  right-left  ambiguity.  Travel  time  can 
then  be  linearly  converted  to  range,  via  mean  sound  speed, 
with  high  accuracy  for  intermediate  ranges  larger  than  a 
few  ocean  depths.  More  complex  conversions  using  depth- 
dependent  slant  range  or  ray  travel  time  are  necessary  for 
local  returns.8  For  synthetic  data,  however,  it  is  possible  to 
maintain  a  spatial  representation  throughout. 

To  generate  synthetic  reverberation  maps,  we  first 
choose  a  scattering  strength  representation  of  the  ocean 
basin  as  defined  by  a  single  arbitrary  image,  shown  in  Fig. 
I.9  We  can  use  a  single  image  because  we  have  assumed 
that  scattering  strength  is  independent  of  incident  and  ob¬ 
servation  angle. 
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The  image  is  composed  of  discrete  pixels  that  form  a 
Cartesian  array  of  grid  points.  We  arbitrarily  choose  scat¬ 
tering  strength  to  be  equal  to  the  depth  value,  in  hectome¬ 
ters,  of  respective  pixels  in  the  image.  This  serves  two  pur¬ 
poses.  It  provides  scattering  strength  values  that  fall  within 
the  range  typically  measured10  and  it  describes  a  scattering 
region  that  has  variations  related  to  the  bathymetry. 

We  define  the  scattering  coefficient  at  each  grid  point 
by  a  discrete  function  a[m9n]9  which  samples  the  continu¬ 
ous  variable  a(x9y)  according  to 

Nx 

cr[m9n]=  X  X  <5(x— m  &x)8(y— n  Ay)<j(xj>). 

m= 1  n=\ 

(1) 

The  scattering  coefficient  a[m9tt\  is  related  to  scattering 
strength  via  2[m,«]  =  10  log(<r[«,m])  for  grid  point  m=x/ 
A x9  n  =y/ Ay,  where  (x,y)  and  [m9n]  represent  respective 
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continuous  and  discrete  (east,  north)  coordinates  originat¬ 
ing  at  the  southwest  comer  of  the  sampled  region. 

An  ideal  line  array’s  azimuthal  resolution  is  most  com¬ 
monly  described  by  its  3-dB  beamwidth  13(6).  For  steering 
angles  from  broadside,  O—tt/2,  to  a  transition  angle  6t  near 
endfire,  0=0, 


P(6)=K 


A  1 
L  sin  0 9 


(2) 


where  K  depends  on  the  taper  function  or  window  used. 
We  define  0,  by  the  value  at  which  the  right-hand  side  of 
Eq.  (2)  reaches  the  endfire  beamwidth,11 

P(0)=2A6KjA7L.  (3) 

We  approximate  0(0)  =0(0)  for  O<0<0,. 

The  beamwidth  of  the  simulated  data  is  determined  by 
the  wavelength  A,  the  array  aperture  L,  and  the  taper  func¬ 
tion.  With  A  =  6.1  m,  L=  160  m,  and  AT=  1 .3  for  a  Ham¬ 
ming  window,  this  translates  to  an  azimuthal  resolution  at 
broadside  of  f3(ir/2)  =2.8°.  The  radial  resolution  A r  is  de¬ 
termined  by  the  mean  speed  of  the  deep  sound  channel  c 
and  the  pulse  duration  for  cw  signals  r.  For  c=  1.5  km/s 
and  r=2s,  Ar=cr/2  =  1.5  km.  The  grid  increments  Ax,  Ay 
are  chosen  to  be  equal  to  the  radial  resolution  of  the  sim¬ 
ulated  data,  Ax=Ay=Ar=l.5  km. 

Reverberation  measured  a  distance  r  from  the  array,  at 
an  angle  0  from  its  axis,  is  comprised  of  returns  weighted 
over  azimuth  and  integrated  within  an  annulus  of  width 
A r.  The  azimuthal  weighting  is  determined  by  the  array’s 
beampattem  in  the  0  direction.  In  our  inversion  and  sim¬ 
ulations,  we  ignore  the  contribution  from  sidelobes  in  the 
beampattem,  and  approximate  the  azimuthal  weighting  by 
unity  within  the  3-dB  beamwidth  0(0)  and  zero  outside. 
The  integration  region  is  an  annular  sector  of  area  r0(0)  Ar 
centered  about  the  point  (r,0).  For  our  purposes,  this  in¬ 
tegration  region  defines  the  scattering  area  for  the  temporal 
return  mapped  to  (r,0).  Contributions  from  a  symmetric 
annular  sector  centered  at  (r,—6)  are  also  added  to  ac¬ 
count  for  the  array’s  right-left  ambiguity.  This  is  displayed 
schematically  in  Fig.  2. 

To  generate  a  synthetic  reverberation  map  for  obser¬ 
vation  s,  data  at  grid  points  within  ambiguous  annular  sec¬ 
tors  are  summed  together  incoherently  via 

rrmax  r& max  V  / 

&[w0,/*0]=  2,  Z  8(x-m  Ax) 

3  rmin  ^min  m=\  ft—  1 


X8{y-nAy)Ts{r,0)d0dr,  (4) 


Ts(r,d)  = 


a(r,0) 


fs(r,d) 


cr(r,—d) 


fs(r,-d) 


,  (5) 


where  £)Jm0,w0]  is  the  reverberation  measured  at  grid 
point  m0=x0/Ax,  «0=>’(/Aj.  (The  selected  grid  points  are 
also  shown  in  Fig.  2.)  Polar  and  Cartesian  coordinates  are 
related  by 


x=r  cos(0+<ps)+jci,  .y=rsin(0+<ps)+.yJ,  (6) 

for  observation  s,  array  axis  orientation  q>s,  and  array  cen¬ 
ter  location  (xs,ys).  The  integration  limits  are 


<y-y„> 


FIG.  2.  Synthetic  reverberation  measured  at  (r,8)  from  an  array  posi¬ 
tioned  at  (x stys).  The  heavy  dotted  line  indicates  the  array  axis  at  angle 
<ps  with  respect  to  the  Cartesian  grid.  Symmetric  annular  sectors  about  the 
array  axis  encompass  grid  points  to  be  integrated. 

emin=0o-p(d)/2,  0mai=0o+0(0)/2,  rmjn=r0—Ar/2, 
^max  roH"  At/2. 

The  function  /s(r,0)  consists  of  two  factors 
fs(r,6)=gs(r,6)hs(r,6).  The  first  factor  gs(r,6)  incorpo¬ 
rates  two-way  propagation  between  source  and  scatterer. 
For  monostatic  situations,  reciprocity  is  used  on  the  return 
trip  so  that  gs(r,6)  =/^(r,0)//sO  where  — 10  log(/5(r,0))  is 
the  transmission  loss  at  the  ocean  bottom  for  horizontal 
position  (r,0)  and  lOlog(I^)  is  the  source  level.  For  in¬ 
version  with  actual  data,  the  insonifying  field  is  computed 
using  a  range-dependent  propagation  model  such  as  the 
wide  angle  parabolic  equation. 

However,  for  this  paper  we  have  assumed  zero  trans¬ 
mission  loss  and  source  level,  i.e.,  gs(r,0)  =  1  for  all  (r,0) 
and  s .  Weighting  the  data  via  a  complicated  propagation 
model  would  be  lost  on  our  simulations  since  the  weights 
would  simply  be  removed  in  the  preprocessing  stages  be¬ 
fore  the  actual  inversion.  (By  similar  reasoning,  we  neglect 
to  convert  the  data  from  range  to  travel  time  back  to  range 
dependence  by  maintaining  a  range  dependence  through¬ 
out.) 

Due  to  the  right-left  ambiguity  of  the  line-array  data 
we  expect  that  Qs[m0,n0]  =  Qs[mf0,n^  when  r0  =  r'Q,  0O 
=  —  0q  for  observation  s .  The  correction  factor, 

hs(r,6)=rP(6)Ar/i]s[m,n\,  guarantees  that  this  is  so  re¬ 
gardless  of  the  grid  size,  where  m—x/Ax,  n  =y/Ay  are 
related  to  (r,0)  by  Eqs.  (6)  for  observation  s.  Here  the 
number  of  grid  points  per  scattering  area,  or  the  integra¬ 
tion  number,  is 

f'-max  ^ 

Vs[m0,n0]  =  I  XI  8(x—m  Ax) 

3  rmin  3  ^min  fn—  1  ft  =  1 

X8(y-n  Ay)d0dr,  (7) 

where  the  summation  limits  are  as  previously  defined.  This 
factor  is  necessary  because  of  approximations  made  in 
mapping  inherently  polar  reverberation  data  onto  the  Car¬ 
tesian  grid  used  to  digitally  define  the  image.  As  the  grid 
size  decreases  with  respect  to  the  radial  resolution  of  the 
array  Ar,  the  correction  factor  hs(r,6)  approaches  unity 
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and  the  resolution  of  the  inversion  increases.  However,  it  is 
not  advantageous  to  choose  a  grid  size  much  smaller  than 
Ar.  As  the  grid  becomes  finer  than  this,  the  information 
content  remains  constant  but  the  inversion  time  prohibi¬ 
tively  increases. 

For  a  given  observation,  we  do  not  compute  reverber¬ 
ation  for  grid  points  whose  scattering  areas  extend  beyond 
boundaries  corresponding  to  the  edges  of  the  image.  We 
also  do  not  compute  reverberation  for  grid  points  whose 
integration  numbers  exceed  a  chosen  threshold  r ]max .  This 
insures  that  redundant  information  from  endfire  beams, 
which  is  too  poor  in  resolution  to  be  useful,  does  not  slow 
down  the  processing.  However,  thresholding  also  limits  the 
data  available  at  long  range,  which  eventually  limits  an 
inversion’s  effectiveness  at  long  range. 


II.  INVERSION 


We  initialize  our  scattering  coefficient  estimate  cr[m,n] 
at  some  value  within  an  allowable  range 
<  a[m,n\  <  <7max  for  all  grid  points  m  =  l,2,3,...,Nx  and 
n= 1,2,3,. ..,Ny.  For  each  observation  s,  we  construct  rep¬ 
lica  data  Rs[m0,n0\  with  our  estimated  scattering  coefficient 
as  in  Eq.  (4),  replacing  Qs[m0yn0]  with  Rs[m0,n0]  and 
<r[m,/t]  with  a[m,n\.  We  then  compute  a  least-squares  cost 
function,  which  sums  the  difference  squared  between  rep¬ 
lica  data  and  true  data  for  each  observation  and  grid  point, 

S  Nx  Ny 

E{&)  =  n  S  (Qs[m,n]-Rs[m,n])2.  (8) 

5=1  m—  1  n=  1 

Minimization  of  this  cost  function  with  respect  to  the 
a[m,n\  leads  to  a  series  of  S  linear  equations  for  each  grid 
point.  This  can  be  seen  by  expressing  the  replica  data  in 
Eq.  (8)  in  terms  of  a[m,n\.  For  notational  economy  we 
convert  [m,n]  dependence  to  an  index  k =  l,2,3,...,iV,  where 
N=NJVy.  We  also  introduce  a  sifting  function  BSKk 
which  defines  the  beamwidth  in  discrete  space 
k=  1,2,3,. ..,JV  for  observation  s .  It  is  unity  inside  the  scat¬ 
tering  areas  about  fixed  grid  point  k  and  its  ambiguous 
image  point  k ',  and  zero  outside  these  regions: 


d  s  N  N  |2 

0  =  X  E  Qs,k  X  <?*/ S,)<Bs,K,k  t 

a(Jj  5=1  k=  1  K=  1 


(9) 


where  y  =  l,2,3,...,W.  This  leads  to  the  linear  equations 


N 

QsJ=  ^  s,kBs,kJ 

K—  1 


(10) 


for  each  observation  s=  1,2,3,. ..,S  at  each  grid  point 
y  =  l,2,3,...,iV. 

When  the  number  of  equations  coupled  to  a  particular 
set  of  grid  points  is  greater  than  or  equal  to  the  number  of 
grid  points  coupled,  the  equations  for  these  grid  points  are 
properly  constrained.  Satisfaction  of  this  criterion  for  a  set 
of  grid  points  depends  upon  the  observation  geometry,  lo¬ 
cal  resolution,  and  integration  threshold.  Satisfaction  of 
this  criterion  for  all  grid  points  is  more  stringent  than  hav¬ 
ing  the  overall  number  of  nonidentical  equations,  which  is 
at  least  S(N/ 2),  exceed  the  overall  number  of  unknowns, 
N.  [We  note  that  the  number  of  nonidentical  equations  is 


at  least  S{N/2)  due  to  the  Cartesian  grid’s  general  lack  of 
symmetry  about  an  arbitrary  array  axis.  If  all  observations 
s  had  array  axes  symmetric  to  the  grid,  i.e., 
q>s=(s—  1)tt/4,  the  number  of  nonidentical  equations 
would  be  limited  to  S(N/ 2),  due  to  right-left  ambiguity.] 

When  properly  constrained,  these  equations  can  be 
solved  for  scattering  coefficient  ak  by  matrix  inversion.12 
Instead,  we  choose  an  iterative  optimization  procedure,  be¬ 
cause  such  procedures  are  extremely  efficient  in  both  for¬ 
mulation  and  solution  of  problems  involving  a  large  num¬ 
ber  of  intricately  related  parameters.  They  are  ideal  for  our 
situation  since  the  images  that  we  need  to  invert  typically 
contain  between  N=  104  to  106  pixels  and  are  the  result  of 
a  highly  convoluted  process. 

Because  the  problem  is  linear,  we  are  guaranteed  to 
obtain  the  correct  solution  via  an  iterative  random-search 
for  local  minima,  if  the  problem  is  properly  constrained. 
This  can  also  be  seen  in  Eq.  (9).  The  function  we  are 
minimizing  is  parabolic  in  the  parameters  to  be  varied. 
This  insures  that  the  minimum  found  in  the  parameter 
search  space  is  a  global  minimum.  Besides  being  extremely 
efficient,  this  method  is  also  convenient  for  another  reason. 
Since  it  is  essentially  the  same  as  using  simulated  annealing 
at  zero  temperature,  it  is  easily  transformed  to  a  simulated 
annealing  algorithm.  If  there  is  uncertainty  in  the  observa¬ 
tion  geometry,  i.e.,  the  ship’s  location  and  array  heading 
are  not  precisely  known,  these  parameters  can  be  included 
in  the  search.  The  problem  then  becomes  nonlinear  and  a 
simulated  annealing  approach  is  necessary  to  escape  local 
minima  in  the  parameter  search  space. 


FIG.  3.  Synthetic  observation  geometry  with  towed-array  course  indi¬ 
cated  by  the  dotted  lines  and  arrows.  The  coordinate  axis  shown  is  at  the 
center  of  the  region  displayed  in  Fig.  1,  i.e.,  at  mc=xc/Ar=81, 
nc=yc/Ar=&\,  where  tick  spacing  is  Ar.  Monostatic  observation  loca¬ 
tions  are  given  by  asterisks.  Array  axis  orientations  are  within  ±  5°  of  the 
tracks  shown.  Specific  array  positions  and  orientations  (x/Ar,y/Ar,cps) 
are  (mcH-4,«c+ 1,5°),  (mc+2,/ic-f3,550),  (mc— l,/ic+4,114°), 

(mc—  3,/*c+2,133°),  respectively,  for  the  5=  1,2, 3, 4  observations. 
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FIG.  4.  Synthetic  reverberation  maps  10  Iog(Qs)  for  the  s=  1,2, 3, 4  observations.  Outside  the  interior  white  lines,  the  resolution  is  poorer  than  one  grid 
increment.  All  data  to  be  used  for  scatter  coefficient  estimation  are  within  the  exterior  white  line,  i.e.,  poor  resolution  data  from  endfire  beams  and  beyond 
a  threshold  range  are  excluded. 


To  implement  the  inversion  numerically,  we  perturb 
the  estimated  scattering  coefficient  bj  at  a  single  grid  point 
j  by  a  random  multiplicative  factor  e  such  that  d'- 
=  6(7 j ,  and  <7min  <  o'j  <  crmax.  (Specifically,  we  choose  e 
—  io(2*-1)  ,  where  0<^<1  is  a  random  variable.  This  is 
analogous  to  the  method  used  in  Ref.  13  for  additive  per¬ 
turbations.  When  ebj  <  amin  or  amax  <  ebj  we  set  d'  equal 
to  Omin/ebj  or  <r^ax/edy ,  respectively.)  The  scatter  coeffi¬ 
cient  bk  remains  the  same  for  all  other  grid  points,  i.e.,  for 
k— 1,2,3. ,.,N  but  k^j.  We  then  compute  a  new  cost  func¬ 
tion  E\  which  is  identical  to  E  in  all  terms  except  those 
containing  the  grid  point  j.  When  k= j ,  bj  is  replaced  by 
the  perturbed  value  d'- .  If  the  perturbed  cost  function  Ef  is 
less  than  the  original  cost  function  E  we  accept  the  pertur¬ 
bation,  i.e.,  if  AE=E'—  E<0.  Otherwise  we  retain  the 
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original  value.  We  continue  this  process  until  each  grid 
point  has  been  perturbed,  i.e.,  for  j  — 1,2,3.. .,N.  That  com¬ 
prises  the  first  global  iteration.  We  then  iterate  globally 
until  the  cost  function  reaches  a  minimum  value  Emin  with 
respect  to  the  bj .  As  discussed  in  the  last  paragraph,  the 
resulting  bj  will  be  the  solution  for  a  properly  constrained 
problem.  If  the  problem  is  not  properly  constrained  for 
certain  grid  points,  resolution  will  be  poorer  than  the  grid 
size  and  speckling  will  occur  in  that  portion  of  the  image. 
Also,  ambiguity  will  be  reduced,  but  may  not  be  elimi¬ 
nated. 

In  comparing  Ef  to  E  it  is  inefficient  and  prohibitively 
time  consuming  to  compute  each  cost  function  as  shown  in 
Eq.  (8).  Since  only  a  single  test  grid  point  has  been 
changed  there  is  no  difference  between  many  correspond- 
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FIG.  7.  Scattering  strength  estimate  2  from  a  dB  average  of  reverberation 
FIG.  5.  The  number  of  reverberation  maps  available  per  grid  point,  fl,.  maps. 


ing  terms  in  each  summation.  Therefore,  instead  of  sum¬ 
ming  over  all  Nx  and  Ny  grid  points  for  each  perturbation, 
only  a  smaller  influence  region  %J[m— m0in— n0]  about  the 
test  grid  point  [m0,«0]  need  be  considered.  In  general,  the 
influence  region  is  approximately  equal  to  the  integration 
region  of  Eq.  (7)  for  a  given  observation  s.  For  observation 
points  related  by  the  array’s  right-left  ambiguity,  the  in¬ 
fluence  region  must  be  large  enough  to  account  for  the  lack 
of  perfect  reciprocity  about  the  array’s  axis.  As  noted  ear¬ 
lier,  this  occurs  for  array  orientations  which  are  not  sym¬ 
metric  with  respect  to  the  Cartesian  grid. 

The  running  time  of  the  algorithm  for  each  global  it¬ 
eration  is  then  proportional  to 

S  Nx  Ny  Nx  Ny 

E  I  I  I  I  &[m-m0,«-Mo].  (11) 

s—  1  /wq=  1  n$=  1  m=  1  n=  1 

This  reduces  to  a  minimum  INflyS  when  the  resolution  is 
equal  to  the  grid  size  for  all  observations  and  grid  points, 
i.e.,  when 


FIG.  6.  Scattering  strength  estimate  2  from  a  linear  average  of  reverber¬ 
ation  maps. 
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%s[m-m0,n-n0]  =8[m-m0in-n0\ 

+S[m-mo(s),«-«o(s)] 

for  all  s9  n0,  m0,  where  [mQ(s),n'0(s)]  is  the  ambiguous 
point  for  [n0>mQ]  and  observation  s.  This  minimum  run 
time  per  global  iteration  is  equal  to  the  overall  run  time  of 
an  average  of  reverberation  maps.  As  we  will  show  in  the 
next  section,  however,  averaging  generally  gives  an  unsat¬ 
isfactory  estimate. 

III.  SIMULATION  RESULTS 

A.  Observation  geometry  and  reverberation  maps 

Simulations  are  preformed  for  the  typical  monostatic 
observation  geometry  sketched  in  Fig.  3.  The  ship’s  course 
is  indicated  by  dashed  lines  and  arrows.  The  straight  line 
segments  form  a  star  pattern,  consisting  of  four  towed- 
array  tracks.  Asterisks  indicate  the  location  of  the  obser¬ 
vations,  where  the  array  axis  is  roughly  parallel  to  the 
track.  To  demonstrate  the  generality  of  the  method,  the 
tracks  are  chosen  to  be  asymmetric  with  respect  to  the 


FIG.  8.  The  cost  function  E,  versus  global  iteration  i  for  the  scattering 
coefficient  inversion. 
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TABLE  I.  Statistical  comparison  between  estimated  scatter  coefficients  a  and  original  coefficients  cr.  Here,  the  estimate  results  from  a  linear  average  of 
reverberation  maps.  Statistics  are  computed  separately  for  each  of  four  decades  in  a  and  for  the  same  respective  grid  points  in  a.  Statistics  are  only 
computed  for  the  21  055  grid  points  where  1  <ft,  see  Fig.  9.  The  percentage  of  this  number  for  each  decade  is  indicated. 


Linear  average 

10-5«7<10“4 

10-',<(7<10-3 

10_3<a<  10-2 

i<r3<<T<io-' 

% 

32.30 

60.81 

6.78 

0.11 

Poo 

0.041 

0.33 

0.88 

0.28 

14.19 

5.01 

2.58 

2.11 

Cartesian  grid.  Also,  to  allow  for  deviations  typical  in  real 
experiments,  the  array  axes  have  been  randomly  perturbed 
from  the  track  by  rotations  of  within  =fc  5°. 

A  synthetic  reverberation  map  10  log (Qs[m,n]),  where 
m=  1,2,3,...,#*,  n  =  1,2,3,...,#^,,  for  each  observation, 
s=  1,2, 3,. ..,5,  is  shown  in  Fig.  4.  Note  the  decrease  in  res¬ 
olution  for  beams  approaching  endfire  and  for  distant  scat¬ 
tering  areas,  as  well  as  the  right-left  ambiguity  about  the 
array’s  axis.  The  interior  white  contour  encompasses  data 
with  spatial  resolution  equal  to  the  grid  increment,  i.e.,  the 
scattering  areas  contain  a  single  grid  point  so  that  77=1. 
Outside  this  contour,  resolution  is  poorer,  i.e.,  the  scatter¬ 
ing  areas  contain  multiple  grid  points  and  1  <  17.  Due  to 
discretization,  these  contours  vary  with  observation  loca¬ 
tion  and  ambiguous  lobes  are  not  symmetric  with  respect 
to  the  array  axes.  The  exterior  white  contours  encompass 
data  to  be  used  in  subsequent  scatter  coefficient  estimation. 
These  data  are  within  the  selected  T7max=  10  integration 
number  threshold  and  are  also  limited  by  a  maximum 
range  from  each  observation  location  such  that 

Rm&x=$3Ar=  124.5  km,  which  spans  about  two  deep 
ocean  convergence  zones.  Data  outside  these  borders  are 
considered  to  have  resolution  too  poor  to  be  useful. 

Due  to  thresholding,  the  mapping  Qs[m0in0]  may  not 
exist  at  grid  point  [m0,n0]  for  observation  s.  To  quantify 
this  for  the  given  observation  geometry,  the  number  of 
mappings  per  grid  point  ft[m,w]  is  shown  in  Fig.  5. 


B.  Estimate  via  linear  average 

We  first  show  the  scattering  strength  estimate  resulting 
from  a  linear  average  of  reverberation  maps.  We  add  all 
mappings  Qs[m0,n0 ]  of  grid  point  [m0,n0]  and  divide  by  the 
number  of  mappings,  0 <ft[m0,/i0]<S,  for  that  grid  point. 
The  resulting  scattering  coefficient  estimates  a[m,n]  are 
converted  to  scattering  strengths  via 
=  10  log {a[m,n])9  for  m  =  1,2,3,...,#*,  and  n=\,2,3,...,Ny. 
The  resulting  image,  in  Fig.  6,  bears  a  poor  resemblance  to 
the  actual  scattering  strength  image  in  Fig.  1.  In  general, 
the  magnitude  and  shape  of  estimated  scattering  features 


are  significantly  different  from  the  actual  ones.  Artificial 
features  of  high  scattering  strength  also  appear  due  to  im¬ 
perfect  removal  of  right-left  ambiguity. 

To  quantify  these  differences,  and  so  test  the  accuracy 
of  the  estimate,  we  compute  the  ratio  of  the  means  and 
normalized  covariance  for  the  original  and  estimated  scat¬ 
tering  coefficients.  We  compute  these  statistics  indepen¬ 
dently  for  sets  of  grid  points  corresponding  to  10-dB  scat¬ 
tering  strength  bins  in  the  original  image,  and  only 
examine  grid  points  where  the  number  of  mappings 
>  1.  The  results  are  in  Table  I,  where  we  use 

<(<?-<£>)  (£7- <CT>)>  ,1<1X 

to  define  the  normalized  covariance  so  that  0<p^a<l.  We 
see  that  the  correlation  is  only  high  for  a  small  portion  of 
the  region  imaged,  where  high  amplitude  scatter  coeffi¬ 
cients  are  found.  The  mean  value  is  estimated  best  for  these 
high  amplitude  coefficients,  but  is  still  only  within  a  factor 
of  3  of  the  true  coefficients. 

C.  Estimate  via  dB  average 

We  next  compute  the  scattering  strength  estimate  re¬ 
sulting  from  a  dB  average  of  reverberation  maps,  i.e.,  we 
add  all  mappings  \0\og(Qs[m,n\)  and  again  divide  by 
il[m,n].  The  resulting  image,  Fig.  7,  has  better  right-left 
ambiguity  reduction  than  the  linear  estimate,  but  still  re¬ 
sembles  the  original  image  only  poorly.  The  comparison  is 
again  made  quantitative  by  computing  cross  statistics,  as 
indicated  in  Table  II.  While  correlations  have  generally 
increased,  they  are  still  low,  except  for  a  small  fraction  of 
the  region  corresponding  to  the  high  amplitude  scatterers. 
However,  the  mean  estimate  is  orders  of  magnitude  worse 
than  in  the  linear  average. 

D.  Estimate  via  global  inversion 

The  global  inversion  method  provides  a  solution  with 
minimum  cost  function  Emin^E4000=EQ(2.17XlO~5)9 
where  the  corresponding  global  iteration  number  i  appears 


TABLE  II.  Statistical  comparison  between  dB  average  estimate  and  original  scatter  coefficients.  See  Table  I  caption  for  further  details. 


dB 

Average 

l<r5<£7<10-4 

l<T4<ff  <  10~3 

10-3<c7<  10-j 

10-3<a<  10_l 

% 

32.30 

60.81 

6.78 

0.11 

Poo 

0.059 

0.53 

0.81 

0.26 

<*>/M 

9.25  XlO-4 

8.34X10"4 

5.27X  10“3 

1.68X10-2 
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as  a  subscript.  Here,  the  cost  function  is  given  in  terms  of 
the  initial  “energy”  E0>  which  is  obtained  by  setting  the 
replica  data  Rs[m,n]  to  zero  in  Eq.  (8).  The  cost  function 
is  plotted  as  a  function  of  global  iteration  in  Fig.  8. 

The  scattering  strength  estimate  for  the  global  inver¬ 
sion  is  shown  in  Fig.  9,  where  the  white  contour  contains 
grid  points  where  1  <  O.  As  expected  from  the  dramatic 
decrease  in  the  cost  function,  the  global  inversion  image  is 
nearly  identical  to  the  true  image  in  Fig.  1.  The  right-left 
ambiguity  of  the  data  has  been  removed  and  the  ocean 
basin  is  resolved  to  the  grid  increment  in  areas  where  re¬ 
verberation  mappings  of  such  high  resolution  do  not  exist. 
This  is  confirmed  quantitatively  in  Table  III,  where  near 
perfect  correlation  is  found  between  the  inversion  estimate 
and  the  actual  scatter  coefficients.  The  mean  values  are 
identical. 

We  note  that  both  high  and  low  values  are  found  with 
equal  precision  in  the  inversion.  However,  low  values  are 
more  statistically  sensitive  to  the  successive  approxima¬ 
tions  since  they  are  given  less  weight  in  the  cost  function, 
when  considered  individually.  As  a  result,  more  iterations 


may  be  necessary  to  reduce  the  variance  of  the  estimate  for 
lower  values.  This  is  shown  in  Fig.  10  and  Table  IV,  where 
the  scatter  coefficient  estimate  and  cross  statistics  are  pro¬ 
vided  for  intermediate  iterations.  For  inversions  with  real 
data,  this  means  that  if  the  cost  function  reaches  an  early 
minimum,  i.e.,  it  is  reduced  by  less  than  an  order  of  mag¬ 
nitude,  a  statistical  analysis  is  necessary  to  interpret  the 
estimate.  (One  possible  way  of  reducing  the  number  of 
iterations  necessary  to  reach  a  minimum  is  to  make  the 
perturbation  decrease  in  magnitude  as  a  function  of  theit- 
eration.  This  has  not  been  done  for  the  present  analysis, 
but  we  have  performed  some  simulations  along  this  line 
which  show  that  the  idea  is  promising. ) 

Even  if  there  are  no  reverberation  mappings  for  a  par¬ 
ticular  grid  point,  due  to  thresholding  or  nearness  to  a 
border,  the  scattering  coefficient  can  still  be  determined  by 
global  inversion  if  the  scattering  areas  of  adjacent  mapped 
points  leak  onto  the  unmapped  point  enough  times.  Com¬ 
paring  the  number  of  mappings  per  grid  point,  Fig.  5,  with 
the  inversion  estimate,  in  Fig.  9,  we  note  that  this  is  the 
case  for  a  variety  of  grid  points  at  the  image  periphery.  It 


dB 


FIG.  9.  Scattering  strength  estimate  2  from  a  global  inversion  of  reverberation  maps,  with  corresponding  global  minimum 
^min~^40oo=^o(2.77xlO-5).  The  white  line  encompasses  grid  points  where  1  <fl. 


990  J.  Acoust.  Soc.  Am.,  Vol.  94,  No.  2,  Pt.  1,  August  1993 


Nicholas  C.  Makris:  Imaging  reverberation  via  inversion  990 


TABLE  III.  Statistical  comparison  between  global  inversion  estimate  and  original  scatter  coefficients.  Estimate  has  corresponding  cost  function 
approximately  at  global  minimum,  Emin = £4000= E0 ( 2.77 X  10"5).  See  Table  I  caption  for  further  details.  (For  grid  points  where  3  <H,p&<7=  1.0  when 
10-5<<7<  10-4.) 


Global 

inversion 

io-5«7<io-4 

10-4«7<10-3 

10-3<ff<10"2 

10-2<<t<  10_i 

% 

32.30 

60.81 

6.78 

0.11 

Poo 

0.93 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

is  also  clear  that  the  leakage  is  insufficient  to  properly  es¬ 
timate  scattering  strength  for  some  other  peripheral  grid 
points. 

IV.  DISCUSSION 

We  have  developed  a  method  for  imaging  the  ocean 
basin  via  a  simultaneous  inversion  of  multiple  reverbera¬ 


tion  measurements  made  at  differing  spatial  locations.  By 
optimal  use  of  the  data,  this  method  removes  the  right-left 
ambiguity  associated  with  line-array  data  and  produces  as 
fine  a  resolution  as  possible.  The  method  has  been  tested 
with  a  typical  monostatic  observation  geometry  where  it 
has  proven  to  be  far  superior  to  data  averaging  for  simu¬ 
lations  involving  reverberation  from  scattering  sites 


dB 

-17.6 


-33.3 


-49.0 


E202=eox10~3  E640=E0x  1CT4 

FIG.  10.  Intermediate  scattering  strength  estimates  2  from  a  global  inversion  of  reverberation  maps  for  indicated  iterations. 
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.  TABLE  IV.  Statistical  comparison  between  four  intermediate  global  inversion  estimates  and  original  scatter  coefficients.  See  Table  I  caption  for  further 
details. 


Intermediate 

iterations 

10-5«r<10-4 

10-4<ct<  10-3 

10-3<cr<10-2 

10~2«7<10-1 

% 

32.30 

60.81 

6.78 

0.11 

(^12  =  -^oX  10  *) 

Paa 

0.030 

0.33 

0.44 

0.17 

(a)/(cr) 

2.09 

1.15 

0.80 

0.49 

(E29=E0X  10-2) 

Paa 

0.12 

0.65 

0.80 

0.39 

(a)/(<r) 

1.26 

1.01 

0.98 

0.86 

(E2(y2=E0X  lO"3) 

Paa 

0.47 

0.95 

1.00 

0.97 

(dr)/ (<j) 

1.02 

1.00 

1.00 

0.99 

(^640=^oXlO"4) 

Pirn 

0.78 

0.99 

1.00 

1.00 

( ar)/(o ) 

1.00 

1.00 

1.00 

1.00 

broadly  distributed  in  space  and  scattering  strength  mag¬ 
nitude.  In  situations  where  only  isolated  scatterers  of  high 
amplitude  are  present,  and  these  do  not  interfere  in  ambig¬ 
uous  beams,  data  averaging  may  suffice.  However,  this  sit¬ 
uation  is  the  exception  for  most  ocean  environments.  A 
more  robust  method,  such  as  the  inversion  presented,  is 
generally  necessary. 

While  the  method  is  currently  applicable  to  many 
long-  and  short-range  monostatic  and  bistatic  experiment 
geometries,  it  cannot  be  applied  arbitrarily.  For  example,  a 
series  of  separate  inversions  is  necessary  to  determine  the 
angular  dependence  of  bottom  reverberation.  In  each  in¬ 
version,  the  separation  between  observations  must  be  small 
enough  to  maintain  a  similar  orientation  to  the  region  be¬ 
ing  imaged.  This  insures  that  the  differences  between  ob¬ 
servation  angle  are  small  with  respect  to  variations  in  the 
angular  dependence  being  measured. 

Alternatively,  it  is  possible  to  reformulate  the  inver¬ 
sion  to  search  directly  for  regional  angular  dependence  via 
a  single  global  inversion.  However,  this  approach  runs  into 
difficulties  for  ocean  environments  dominated  by  negative 
excess  depth  and  rugged  bathymetry  due  to  extreme  vari¬ 
ations  in  transmission  loss  that  can  make  spatially  diver¬ 
gent  observations  independent  of  each  other.4  For  example, 
regions  of  extremely  high  transmission  loss,  or  shadow 
zones,  may  exist  where  no  useful  scattering  information 
can  be  extracted.  Under  these  circumstances  scaling  by 
transmission  loss  is  not  appropriate  and  redundant  infor¬ 
mation  vital  to  ambiguity  resolution  may  be  lost.  A  simple 
way  to  circumvent  this  problem  is  to  invert  observations 
sufficiently  close  that  their  respective  transmission  loss 
functions  are  approximately  the  same  and  to  exclude  data 
when  transmission  loss  exceeds  a  threshold  maximum.  Op¬ 
timally  resolved  linear  reverberation  would  then  be  the 
output  of  the  inversion  rather  than  scattering  coefficient. 

On  the  other  hand,  variations  in  transmission  loss 
along  ambiguous  radials  can  be  used  more  advantageously 
than  has  been  suggested.  These  variations  often  provide 
additional  environmental  symmetry  breaking  information 
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that  can  be  exploited  to  help  resolve  right-left  ambiguity 
for  a  single  observation,  as  is  demonstrated  in  an  upcoming 
article.4  The  present  inverse  method  can  be  readily  modi¬ 
fied  to  exploit  these  situations,8  and  provide  a  means  of 
directly  estimating  the  regional  angular  dependence  of 
scattering  via  a  single  global  inversion  of  spatially  diver¬ 
gent  observations. 

No  attempt  has  been  made  to  account  for  statistical 
variations  in  reverberation  that  may  arise  at  differing  ob¬ 
servation  locations.  While  such  stochastic  phenomenon  are 
possible  in  real  data,  they  are  not  relevant  to  this  initial 
formulation.  Our  goal  is  principally  to  show  that  there  are 
ways  to  make  optimal  use  of  towed-array  data  that  are 
independent  of  the  often  stochastic  nature  of  the  reverber¬ 
ation.  And,  even  if  high  variance  exists,  an  ensemble  of 
observations  with  redundant  geometry  may  be  used  in  the 
present  formulation  to  make  the  problem  statistically  over¬ 
determined  and  so  reduce  the  variance  by  the  number  of 
redundant  observations.  Alternatively,  stochastics  can  be 
incorporated  by  searching  for  probability  distributions  to 
parametrize  seafloor  scattering  rather  than  mean  scattering 
coefficient  values.  If  the  reverberation  is  inherently  deter¬ 
ministic,  i.e.,  the  scattering  surfaces  are  much  larger  than  a 
wavelength  and  approximately  planar,  or  the  variance  is 
sufficiently  low,  only  observations  sufficient  to  remove  am¬ 
biguity  and  obtain  the  desired  resolution  would  be  re¬ 
quired,  as  in  the  simulations  presented. 

Finally,  we  contend  that  it  is  possible  to  determine  the 
optimal  set  of  observations  needed  to  image  the  ocean  bot¬ 
tom  with  methods  related  to  those  presented.  That  is,  as¬ 
suming  some  or  no  knowledge  of  regional  bottom  scatter¬ 
ing,  and  given  constraints  upon  our  resolution,  ship 
maneuverability,  time,  site  locations,  etc.,  we  could  invert 
for  the  optimal  course  of  the  research  vessel.  This  has  been 
referred  to  as  the  “traveling  acoustician  problem.”  14  Qual¬ 
itatively,  the  problem  is  related  to  the  much  simpler  prob¬ 
lem  of  refining  the  array’s  navigation  by  including  it  as  a 
free  parameter  in  the  reverberation  inversion.  Both  prob- 
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lems  are  nonlinear,  and  are  most  efficiently  solved  with  a 
technique  like  simulated  annealing. 
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